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Abstract 

We consider a new kinetic equation for systems with a multistep potential of inter- 
action proposed by us recently in Physica A 234 (1996) 89. This potential consists 
of the hard sphere part and a system of attractive and repulsive walls. Such a 
model is a generalization of many previous semi-phenomenological kinetic theories 
of dense gases and liquids. In this article a normal solution to the new kinetic equa- 
tion has been obtained, integral conservation laws in the first order on gradients of 
hydrodynamic parameters have been derived as well. The expressions for transport 
coefficients are calculated for the case of stationary process. We also consider lim- 
iting cases for this kinetic equation. For specific parameters of model interaction 
potential in shape of the multistep function, the obtained results rearrange to those 
of previous kinetic theories by means of the standard Chapman-Enskog method. 
In view of this, new theory can be considered as a generalized one which in some 
specific cases arrives at the results of previous ones and in such a way displays 
the connection of these theories between themselves. At the end of this article we 
present results of numerical computation of transport coefficients for Argon along 
curve of saturation and their comparison with experimental data available and MD 
simulations. 
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1 Introduction 

The main problem of kinetic theory in the last few decades was to appreciate 
properties of dense gases and liquids on the basis of motion laws and molecular 
interactions. In 1872 L.Boltzmann suggested a kinetic equation for a dilute gas 
[1]. At the same time, he constructed the entropy functional and proved H- 
theorem. 

In 1946, N.N.Bogolubov suggested a new approach [2] for the problem of 
derivation of kinetic equations, which describes irreversible processes in gases. 
He aspired to change suggestions of the Boltzmann's theory to a more rigorous 
ones and generalize the theory to higher densities. In this theory the Boltz- 
mann kinetic equation is obtained in the zeroth approximation on density. It 
is possible at the description of gas kinetics on rough enough time scale and 
at the presumption that quantities, which define gas state, can be expanded 
into a series on the density. However, after studying of higher order expansion 
terms (starting from the second one) their divergency has been found [3] . The 
reason for such a divergency lies in the fact that the accepted expansion is 
based on a dynamics of isolated groups of particles in an infinite space with- 
out taking into consideration the whole media, i.e. without taking into account 
all others particles. In contrast to an equilibrium liquid, velocity correlations 
in nonequilibrium state include long-range correlations between particles. It 
does not allow to generalize the Boltzmann kinetic equation to high densities. 

It is particularly remarkable that a successful empirical kinetic theory of dense 
gases was developed by Enskog in 1922 at least for the case of hard spheres 
(standard Enskog theory - SET) [4]. His arguments were similar to the Boltz- 
mann's ones. In 1961 Davis et all [5] suggested the DRS kinetic theory, where 
the so-called "square-well" potential of interaction is considered. DRS theory 
is the analogue of the usual SET theory for a specific type of interaction poten- 
tial. Here, an attractive part of real interaction potential is approximated by 
some finite attractive wall. A revised Enskog theory RET [6-8] and a revised 
version of DRS - RDRS [9] have been obtained. The necessity of the revised 
versions could be explained as follows. I: kinetic equations of initial versions of 
those theories cannot be derived sequentially. II: exact entropy functional can- 
not be constructed, as a result //-theorem cannot be proved. But some recipe 
of construction of entropy functional for the SET theory was suggested in pa- 
pers [10-12] The exact entropy functional and proof of the i/-theorem were 
constructed for revised versions RET [7,8] and RDRS [9] only. Furthermore, 
the RET kinetic equation has been obtained successfully within the frame of 
some theoretical scheme, which is analogous to Bogolubov one [2], but using 
a modification of boundary conditions. The last ones take into account local 
conservation laws in the solution of the BBGKY hierarchy [13]. 

For the purpose to use results of the SET theory to real systems, Enskog 
suggested to change hydrostatic pressure of a system of hard spheres to a 
thermodynamic pressure of a real system. Having this assertion in mind, Han- 



ley et all [14] built new kinetic theory MET (modified Enskog theory), where 
hard sphere diameter a is defined via the second virial coefficient of the equa- 
tion of state of a system. In such a way, a becomes dependent on temperature 
and density. Using different equations of state: BH [15], WCA [16], MC/RS 
[17,18] and others, one obtains corresponding versions of MET. 
In the kinetic mean ffeld theory (KMFT) [19] next to the hard sphere interac- 
tion potential one considers also some smooth attractive "tail" . It is noted in 
[20] that in this case the quasiequilibrium binary correlation function of hard 
spheres should be replaced to a such one, which takes into account explicitly 
or implicitly the total interaction potential. The main conclusion of KMFT 
is that the smooth part of interaction potential in the first approximation 
on gradients of hydrodynamic parameters does not contribute explicitly into 
transport coefficients. There is only an indirect contribution via the binary 
correlation function. Its dependency on temperature is defined by a smooth 
part of the interaction potential. 

In our recent paper [21] we suggested a new kinetic equation for systems with 
a multistep potential of interaction (MSPI). This potential consists of the hard 
sphere part and a system of attractive and repulsive walls. Such a model is 
a generalization of SET (RET, MET), DRS (RDRS) and KMFT theories. 
We also proved for this equation if-theorem. But normal solution was not 
published yet. In this article going to the schema of construction of normal 
solutions of kinetic equations with the help of boundary conditions method 
[22] , a normal solution to the new kinetic equation has been obtained, integral 
conservation laws in the first order on gradients of hydrodynamic parameters 
have been derived as well. The expressions for such transport coefficients as 
bulk and shear viscosity and thermal conductivity are calculated for the case 
of stationary process. We also consider limiting cases for this kinetic equation. 
For specific parameters of model interaction potential in shape of the multistep 
function, the obtained results rearrange to those of the SET (RET, MET), 
DRS (RDRS) or KMFT theories by means of the standard Chapman-Enskog 
method [23]. In view of this, new theory can be considered as a generalized 
one which in some specific cases arrives at the results of previous ones and in 
such a way displays the connection of these theories between themselves. In 
section 9 we present results of numerical computation of transport coefficients 
for Argon and their comparison with experimental data available and MD 
simulations. 

2 Multistep potential of interaction. Kinetic equation 

Let us consider a system of A^ classical particles in volume V when N -^ oo 
and \/ — » oo provided N/V = const. The Hamiltonian of this system reads: 

N 2 

HN = Y.^ + Y.V^r (2.1) 



Our purpose is the most detailed analysis of nonequilibrium processes in dense 
systems. To do this let us define MSPI ipij = ip{\ri — rj\) = ip{\fij\) = ^{rij) 
in a form of the multistep function: 
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Here A^* is the total number of attractive and repulsive walls except hard 
sphere one. For our convenience we distinguish separately systems of attractive 
and repulsive walls. Let one has n* attractive walls, they are separated by the 
distances an and have heights Aeu, i = 1, . . . , n*; and m* repulsive walls with 
parameters arj and Ae^j, j = l,...,m*, correspondingly. o"o is the location 
of hard sphere wall. It is obviously that n* + m* = N*, Aeu = ^h — £ii+i, 
ASrj = Erj+i — Srj. lu such a Way, parameters ao, n*, an, Aeu, ^, CTrj, ASrj 
define multistep potential of interaction completely (see Fig. 1). 




Fig. 1. Multistep potential of interactin with parameters re* = 1, m* = 3. 

Going similarly to the derivation of the kinetic equation of the RET theory 
[7,8,24] and taking into account the system of attractive and repulsive walls, 
one obtains the following kinetic equation: 



( d d \ f - 

I at ^'"^Wl •^i^^i'^^^y dx2T/2(xi,X2;t), (2.3) 

/2(xi, X2] t) = gj{fi, f2|n(t), /5(t))/i(xi; t)/i(x2; t). 



where g2(ri,r2\'n'{t), P{t)] is defined in the usual way from maximum of the 
entropy functional and in its turn is the functional of local values of density 
n{fi; t) and temperature I3{ri\ t). T is an operator which describes interaction 
of two particles at presence of MSPI: 



f = f^, + Y.Tu + Y.Tr,^ 



(2.4) 



i=l 



i=i 



fj^ = al J da ag d{ag) {5 [n - ?2 + -^^0+) fi„(a) - 5 (n - r^ - 0-^0+) } . 

(2.5) 

The last expression is nothing but operator of hard spheres interaction [25], 
a is the unit vector along distance from the second particle to the first one, 
g = V2 — v\ Ys, relative velocity. Baia) is the velocities shift operator as in 
classical mechanics of elastic collisions. Tu = J2k=b,c,dTiik is an interaction 
operator on ith repulsive wall; Trj = J2k=b,c,dTrjk is an interaction operator 
on jth attractive wall; a, b, c are types of possible interactions (see Fig. 2). 
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Fig. 2. Types of possible interactions (schematic draw). 
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Que = e(ag- —^j , Or.c = e ^-ag - — -^ j , (2.7) 

Ou, = 9 f i^ - ag) 9iag), 9,.,, ^ 9 f i^ + ag) 9{-ag). 

Expressions (2.7) are nothing but conditions for specific type of interaction 
{9{z) is the Heaviside unit step function). 

B " ^(o") is an operator which acts on velocities and changes them in accordance 

with the interaction of specific type on each waU: 

Bu^{a)<iJ{v,,V2) = <il{v^^„v^i^). (2.8) 
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type "d": 

<'' =vi + cj{ag)=v',. (2.12) 
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Now kinetic equation (2.3) can be rewritten in an exphcit form: 

I + ajr) /■(...)= (2.13) 
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JeUi, fl) + 2^ (Jlib + Jlic + </iid) + 2_^ {Jrjb + Jrjc + Jrjd)- 
i=l j=l 

One can draw the following conclusions regarding the form of (2.13). In the 
absence of attractive and repulsive walls (Aeu = 0, i = 1, . . . ,n*, ASrj = 0, 
j = 1, . . . , m*) this kinetic equation transfers to that one of the RET theory 
[7,8]. In the presence of only one finite attractive wall (Aeu = 0, i = 1, . . . , n*, 
AErj = 0, j = 2, . . . ,m*, AEri 7^ 0) one obtains the kinetic equation of the 
RDRS theory [9]. Moreover, it can be shown, in the third special case when the 
set of walls is merged with some smooth potential (pt andACT^j = cru — au^i -^ 0, 
i = 1, . . . ,n* — l,Aeii — > 0, n* ^ oo, Adrj = (Trj+i — cTrj ^ 0, j = 1, . . . , m* — 1, 
AErj ^ 0, m* -^ oo, and 
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kinetic equation (2.13) transfers to those one of the KMFT theory [19]. 



3 Macroscopic conservation laws. The zeroth approximation 

Let us introduce the following set of hydrodynamic parameters: particles num- 
ber density n{r;t), hydrodynamic velocity u{r;t), densities of kinetic Wk and 
interaction ooi energies: 



n(fi;t) = / dt/s fi{xi;t), 

f Vl 

u{fi]t) = / dt/s fi{xi]t) 



(3.1) 
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o.k(r,;t)=/d.lA(a:,;t)^^, 

uj;{fi-t) = Jdf2n{f2;t)g^{fi,f2\n{t),P{t))—<^{\fu\). 

We also introduce the set of interaction invariants ^=|m, mv, ^frnf ^+0(ri2)]k 

Multiplying initial kinetic equation (2.13) on each component of the vector \E' 
and integrating with respect to V2 one obtains the equation of continuity, the 
equation of motion and equation of balance of kinetic energy, correspondingly: 
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dt mn [ or la or 113 J 

The equation for interaction energy density one obtains automatically differ- 
entiating the expression for uji taking into account conservation laws for n{r\ t), 
u{f;t) and uy^ (3.2). Stress tensor P^js and heat flow vector consist on kinetic 
and interaction parts: 

PafS = Pi, + PI., , ^33^ 

Qa = QI + Qa ■ 

In its turn, the potential parts consist of the hard sphere term and components 
conditioned by the set of attractive and repulsive walls: 
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(3.4) 



where indices 1, 2, 3, and 4 correspond to different types of interactions: "a", 
"6" , "c" , and "rf" , correspondingly. Now let us write down the expressions for 
stress tensor P^^ and heat flow vector q^ in an explicit form: 



Pap = J dt^i fiixi;t)mciaCifs, 
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Expressions for Pj2ai3, PjSa/s, Pj4ai3, qj2a, qj3a, and qj4a look similar to (3.6)- 
(3.13) Pi2af3, Pi3af3, Pi4af3, qi2a, QiSa, and Qij^a, Correspondingly, at formal re- 
placing a^ -^ a^j, erg -^ -ag, Vn -> Vir, AEu -^ AErj. 

At the end of this section we consider macroscopic conservation laws in the 
zeroth approximation. Let us suppose that one-particle distribution function 
in this case is equal to the local-equilibrium Maxwell one: 



Neglecting by any spatial gradients one obtains 



(3.15) 
1 1 



d\ f2ix,,X2;t) c:. J d\ fn---) = gTia\n,(3)fr{x,;t)frin,V2;t). 



It is well known from the equilibrium statistical mechanics [26] that any sharp 
jump of interaction potential results in a corresponding jump of binary corre- 
lation function: 

cq/ -I ^ - exp IP AElif , 

(3.16) 



g'2^{(^rj\n, (3) 



exp {P AErj} 



Substituting these relations into (3.5)-(3.13), one obtains: Pa/3 = pSa/3, where 
p is the hydrostatic pressure which consists on kinetic p^ and interaction p^ 
parts; qa = 0: 
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p = p^ + p\ 

p^ = n/ceT, 

P' = -Trn^kBTA, ^^^^^. 



A = alg2{a^) - j; alg2{a+) {e-^^^- -l} + f] 4,g2{a;^) {e'^^^^^ -l}, 
g?^g?{a\n{n;t),P{n;t)). 
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Then, starting from (3.2), one has the conservation laws in the zeroth approx- 
imation (Euler laws): 

1 dn du 
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ujf\ri]t) = —Jdfun{fi]t)gl'^(^\fu\-,n{fi,t),f3{fi-,t)^^{\fu\). 



4 The principle of construction of higher order approximations. Normal 
solutions by means of boundary conditions method 

Let us consider initial kinetic equation (2.3) (or (2.13)) 

(^^+v,^] fi{x,;t) = J dx2fMx,,X2;t) = J{h,h), (4.1) 

where the collision integral J(/i,/i) consists of the usual one of the RET 
theory - Je(/i,/i) and collision integrals demanding on interaction on each 
wall. Such a structure of J(/i, /i) is caused by the structure of T-operator: 



J{fl, /l) — JeHi, /l) + Z^[Jlib + Jlic + Jlid] + Z^[Jrjb + Jrjc + Jrjd]- (4.2) 

i=l i=l 

Going to the schema of construction of normal solutions to kinetic equations 
with the help of boundary conditions [22], let us introduce in the right-hand 
side of (4.1) an infinitely small source — ^(/i — /j ), where e -^ 0: 

For deviation 5f = fi — /} this equation reads: 
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J (/;°'. /,'"') + J (/r. if) + J {sf. /f ') + J (sf. if) . 

First and foremost it should be noted that in the colhsion integral in the 
usual Boltzmann equation is local (/i is a function of the same Cartesian 
coordinate ri). In our case the J(/i,/i) is a nonlocal collision integral. One 
time /i is calculated in ri, second time - in ri ± ao", integration with respect 
to a (surface of unit sphere) is performed. As a consequence of that in one 
or another way we will find solutions in some approximation and there is no 
sense "to draw" the whole nonlocal collision integral. It is much convenient to 
use its approximate expression. This approximation should be the same order 
as within the frame for the solution of total kinetic equation. In such a way, 
let us expand functions /i and g2 in the vicinity of ri into series on deviations 
±(Tcr. This results in: 

oo oo 

j = Y.Jk = Jo + J\ r = Y.Jk, (4-5) 

fc=0 fc=l 

oo oo 

/ = ^4 = /o + r, r = Y.h, (4-6) 

fc=0 k=l 

where k is the expansion order. All Jk are here local functionals of /i. We also 
use the following notations: I{5f) = J (fi ,Sf) + J (Sf,fi ) - linearized 
nonlocal collision operator. We have similar expansion for this operator - the 
relation (4.6). Here loi^f) = Jq (/} , 5f) + Jq (5f, /j j is hnearized local col- 
lision operator which coincides with those one of the usual Boltzmann kinetic 
equation within a factor of g2^{crQ) if one neglects the set of walls except hard 
sphere wall only. Then equation (4.4) transfers to 

pf + e6f-I,{6f)= (4.7) 

-^/f ^ + J {fi'\ /f + r{5f) + J{5f, 5f) - v.^^Sf. 

To solve this equation by means the boundary conditions method we need 
in its integral form. To this end let us introduce operator S{t, t') with the 
following properties: 

^^S{tX) = hS{tX). S{t,t) = l. (4.8) 

Using the limiting condition lim 5f{t) = one obtains: 

£— >— oo 



5f{t)= j dt' e'^''-'^ S{tX) X (4.9) 
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Equation (4.9) is completely ready for iteration procedure. This procedure can 
be organized as follows: 

SI<'^'\t) = I dt' e^-(*'-*) S{t, f) l^jf + J^'^+i) {ff\ If) + (4.10) 



d 



r(fe+i) (5/(fc)) + j(fe+i) (<5/(fe)^<5/W) _^^^<5/(fc) 



where 

fc+1 
fc'=0 

Each (/c + l)th step uses conservation laws in kth approximation. 

5 Distribution function /i in the first approximation 

The expression for distribution function /i in the first approximation is ob- 
tained if one puts into (3.5) k = and take into account the equality 5f^'^^ = 0. 
Then we have: 



5/(1) (t)= [ dt' e^(*'-*)5(t,t') 



^/r + ^o(/r,/r) + Ji(/r,/r 



,(0) .(0) 



r-(O) .(0) 



It can be shown that 
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(5.1) 



(5.2) 



Making the expansion with an accuracy of first order on gradients of hydro- 
dynamic parameters 



g^,{n,r,\n{t),P{t))c^gl^ ru-,nr-^-t],P^^-,t 
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/i(fi ± aa) = /i(fi) ± — /i aa 
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and taking into consideration conservation laws (3.18) one obtains after very 
unwieldy calculations the following: 
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4 + — 5t;^ - 9s^ + ^ U^ - s 4t;^ + 6s^ U + 



3/2 



(5. 
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^ e-^^ I dx2 exp {-h' - v'} U-'-^ + l {v' + s'f] 5^p - 

- / (1x2 exp \ —h'^ — f ^ W f ^ 1 — 

9 J •- J \ V V 



v' + s'f"\5. 



Here 



,1/2 



{13 Asii 



{f3Aer,Y^\ 



X2 



m 



4:knT 



1/2 



C2, 



(5.9) 



m 



1/2 



& 



m \V2 ^ - 1 



g =C2-Ci, 
1 

2' 



Then 



(1) _ /(o) 



/i = /r^ = /r^ + sf 



c(l) 



t r 8 



(5.10) 



(5.11) 



d 

■Pa "c 



It can be shown that J dvi 5f'^^\fi, vi, t)^ = 0, i.e. the hydrodynamic param- 
eters n{fi;t), u{fi;t), j3, Ui are completely defined by the local one-particle 
distribution function /} (3.14). 



6 Conservation laws in the first approximation. Stationary process 

First, let us calculate kinetic parts of the stress tensor and heat fiux vector. 
Substituting one-particle distribution function /i (5.11) into (3.5) one obtains: 



Pa^^ =/'^«/3 + J dt' e^(*'-*) M\t, t') [S^p\, 

d d 2 9 



gk(i)= Ut' e<''-'H\tX) 



3(9ri^ 
d 
dru 



InT 



(6.1) 



(6.2) 



Jt' 



where cores of kinetic part of transport laws read: 
M'^it, t') = — J dvi mciaCifsSit, t') {L,^},, , 
L\t, t')=\j dv, c,J^S{t, t') {K^}, . 



(6.3) 
(6.4) 
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For calculation of potential (interaction) parts oi P^^ and q^^^\ the expression 
Jo dA/2(ri + Xaa, Vi,fi + Xaa — aa, V2] t) = Zi + Z2 should be expanded into 
series on, first of all, inhomogeneity of distribution function, then on deviation 
of /} on 5f^^\ In both cases one should keep only the first order on gradients. 
Calculations give: 

^i = 7^f^& W/i '{xi;t)fl >{ri,V2]t)aj—\n ^ (6.5) 

Z an fl'{ri,V2;t) 

Z2 = g?{cT) {/f )(xi; t)<5/(i)(fi, V2; t) + 5/(1) (xi; t)f['\n,V2; t)} . (6.6) 

Zi is the expansion of /| on inhomogeneity (the inhomogeneity of 5f'^^^ is 
considered as a negligibly small quantity), Z2 is the expansion on deviation 
5f^^\ Then, general expressions (3.6)-(3.13) taking into account (6.5) and 
(6.6) transfers to 

i^), = /dt'e^(*'-*)MXt,t')[5.,],- (6.7) 

4 , I r^ „ f6^ d 

1/3 + -^ '"7Wa/3 



n^^JirmkBT H2 I ^S^f^ + -^^u^^af^ } , 



q^ = Jdt' e<''-'^L\t,t') 



^ InT 



<9riQ, 



-n A:B^^^i/.^, (6.^ 



where M\ U are cores of potential parts of transfer laws. Their structure is 
very complicated. To save space, these expressions are not presented here in 
their explicit form. 

n* 

H2 = 4s^{<rt) + E ^uS2{at) e-^^^" S(/? Ae,) (6.9) 

+ E<-&''K-)e-^"^-2(;9A5.,), 



S(s) =e^-\s-Kl{s), 

Kr(s)=2 dxxV 
7o 



x^ + s. 



(6.10) 



In such a way, transport laws in the first approximation are not completely 
integral as this takes place in the case of solution of the usual Boltzmann 
kinetic equation [22]. There are also local-time terms (second terms in (6.7) 
and (6.8)). They are caused by the inhomogeneity of /} and, consequently, 
are not sensitive to the "memory" effects. 

In stationary case operator Iq does not depend on time exphcitly: 



16 



5(t,t')=e^o(*-*'). 
Then 

5/« = y'dr e-e-^«4ir„^ In T + L„^^M J, T = t'-t. (6.11) 

Let us define the following quantities: 

aa{T) =e-''UK^j^\nT\, a„(0) =i^«-^lnT, 

_ f 9 1 (9 ^ ^ ^ 

bap{r) = e ^0^ <^ L„^- — u^ \ , &a/3(0) = L^/j- — ««• 

It can be shown that the operator Jo has the same mathematical properties 
that the corresponding operator of the usual Boltzmann kinetic equation: 

Jo (/r^) = 0, /o(0 = A^, A < 0. (6.13) 

Then 

\\aa{r)\\ ^ ||a«(0)|| exp{K^^T}, 

r < 0, Amax = max|A| (6.14) 

||&a/3(T")|| ^ ||&a/3(0)|| expjAmaxr}, 

and 



V'- 



Using last transformations, equation (6.11) transfers to: 



Sf^'^ = ^hm J dr e- [a^{r) + 6„^(r)) = / dr [a^{r) + 6„^(r)), (6.15) 

— oo — oo 

or, introducing Sf^-^^ = ^'■^-'/i ) it can be rewritten in the following final form: 
<t>^^^ = A^^—\nT + B^p^—u^. (6.16) 

Act and -Bq,/3 satisfy integral equations like this: 

(6.17) 

loiBajs) = Lap- 

Operator Iq has the following structure: 
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Io = Ia + J2 {^^b + I^c + ^id} + E i^Jb + Ijc + Ijd} , (6. 18) 



where 



4(0i) =a2g^<^(ao+) Idt/adS- age{ag) x (6.19) 

J,,(0i) = algl'^iaf^) J dt/sda <Tg e{-ag) X (6.20) 

/f V^l)/f°V"'2) {01 + 02 - 0';. - 0'2'.} , 

hc{<Pi)=crlgTia,:,)ldv2daag9 f^g-^i^ j x (6.21) 

/r(^i)/f°^(^-'2){0i + 02-0';:-0'2':}, 



U<P,)=algl\a^,)Jdv,da ag9{ag)9 U -^ - agj x (6.22) 

/fH^i)/fH^2){0i + 02-0;T-0'2T}, 
J,,(0i) = <g2''«) / d^2da ag ^(ag) X (6.23) 

/f^(^i)/f^(^2){0i + 02-0r,-0'4-}, 
/,,(0i) =<g2''K) /dt/2da ag^ (-ag- ^l^^j X (6.24) 

/f)(t;i)/f°)(tT2){0i + 02-0r;-0^;}, 



lAAEr 



I,Mi)=<^%S^K) j ^vAa age{-ag)e U_^ + ^g] X (6.25) 

f^\v,)f^^\v,) {01 + 02 - 0r; - 0'27} > 



01 = 0(fi,t;i;t), 0^, = 0('fi, {?!',; tj 

02 = 0(fi,tT2;t), 02^ = (fi,t^,;t 



// // /// //// 

Kill 



(6.26) 



In the case of SET (RET) theory, the hnearized local integral operator Iq = la, 
whereas in the case of the Boltzmann kinetic equation there is usual Boltz- 
mann's linearized operator and la tends to that one in the limit of low density: 
n^O,gl\a^)^l,Ia-^lB. 

In such a way, to find one-particle distribution function in the first approx- 
imation in stationary case one should analyse integral equations (6.17) and 
solve them. 



7 Solutions to the integral equations 

To find quantities A^ and Bap we liave set of integral equations (6.17). Let 

US define dimensionless self velocity: w = [2^^] c. Using the property of 
isotropy (in the velocity space) of the operator Jq and structures of Ka (5.7) 
and Lai3 (5.8) solutions to (6.17) can be presented as follows: 

Aa=WiaA{wi), (7.1) 

Ba/3 = Biap + B2ap{Wl)5ap-, (7.2) 

The structure of Ba/s is caused by the structure of Lais'- Lap = Liap + -^2a/3, 
where in Lia/s are all terms with (wiaWii^ — ^wlda/s), in L2af3 are all terms 
with Sai3- Then 

Io{WiaA{Wi)) = Ka, 

Jo iiwiaWi/s - -wlSa^j Bi{wi)] = Li„^, (7.3) 

Io{B2{wi)) = L2. 

Following the standard Chapman-Enskog method [23], let us represent A{wi), 
Bi{wi) and B2{wi) via the Sonine-Laguerre polynomials 



m 



r(n + 171 + 1] 



J^o j!("^-j)'r(n + j + l) 

i.e in the form 



A{w,) =E«^™^^SK 



3/2 

m=0 

oo 

(m) Q{m) ^^^^2 
5/2 

m=0 

oo 



B,iw,)=j:b^r^Sij^iwl), (7.5) 



m=0 

Fredgolm condition gives limitations on expansion coefficients a^'^^ and &("*): 
the correction in the first approximation does not contribute to hydrodynamic 
parameters J dwi /{ (p^-^^ = 0. By this means: 

6(°) = 0, (7.6) 

6« = 0. 
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It is known that the Sonine-Laguerre polynomials converges quickly. Therefore 
only the first nonzero term is considered in expansion. This is as a rule and we 
will follow the procedure. Such an approximation gives the error for transport 
coefficients (practically for all types of interaction potentials) not exceeding 
2%. Thus, we have: 



'5/2 



5 2 



B,{w,) ^ b^l'^SJ^hwl) = 6f), (7.7) 



B,M^b^^^S^y,iw!)=b?^[^-lwl + \wt) 



and 



Jo [wiaa'-'^^ ^ 2 ^^v) ^ "^"' 

Jo ( ( WiaWiis - -wl6at3] &l j = -^1q/3, (7.^ 

,(2) ^15 5„„2 , 1„„4 



Multiplying these equations on wiaSy2i'^i)i ('^ia'y^i/3 ~ I'^i^a/?)) 'S'|,2('"^i); 
correspondingly, and integrating with respect to Wi, one finds: 

(1) ^ J dwi KgWiaS^i/^jwl) 



(o)_ / dwi Li^p (wi^wip - |wf ^a/3) 

1 — 7 i \ — 7 i \"' I ' ^ ^J 

/ &Wi Jo [wiaWip - ^wlSap) [wiaWlfi - 3 t^i^a^ j 

First-hand calculations for a^^'^ and 6*^*') give the following: 



8v^n2 {a* + §^2} 






8v^n2 {a* + ^1^2} V2A;BTy 



where 
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i=l 
m* 

n* m* 

i=i i=i 



ASi = /5Aezi, i = l,...,n*, 
ASj = P AErj, j = I, . . . ,m*, 



A* = - 
2 



i/,(^) = v^/2 + e^r(3/2,z), (7.14) 



r(r, s) = J^ dx x^ ^e ^ is incomplete r-function. As far as 62 does not 

contribute into transport coefficients, we did not calculate it. 

In this manner, one-particle distribution function in stationary case reads: 

/« = /(°)(l + 0«), (7.15) 



where a'-^-* and 6'-^-' are defined by (7.12) and (7.13), correspondingly. 



8 Calculation of transport coefficients. Some limiting cases 

General expressions for the stress tensor and heat flux vector for nonstationary 
process were obtained in section 6. But explicit calculation was performed for 
only one part which is inhomogeneous on /} and local-time next to other 
integral terms. In stationary case integral terms transfers to local-time. The 
explicit calculation for these terms becomes possible. Substituting one-particle 
distribution function /i (7.15) into general expressions (3.5)-(3.13) and taking 
into account new structure for (6.6): 

Z2 = gl''{^)f['\xi;t)f['\n,v,;t) {<P^'\x,;t) + ^^'\n,v,;t)} , (8.1) 

one obtains: 
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d 

PaP=P^aP - l^-Q U^^aH " 2r7S'„;3, (8.2) 

ga = -AT^T, (8.3) 

where S'q,/3 is the velocities shift tensor. Exphcit expressions for transport co- 
efficients bulk K and shear ri viscosities and thermal conductivity A read: 



K = ^n^7rmkBTf^H2, (8.4) 



3 

^=5 



K + ^nkBT h + ^V^nHi\ b^'>\ (8.5) 



3kB ^ ^^'^ ^^^^^ 



2m 
where 



^+lnkJ^] {l + ^v^niJi}a«, 



^.6) 



6(°) = -6f\ ifi = ^A- (8.7) 

E ^f.g?(^.l) A., e-- I ^ - ^ + ^ e-^ r (^, AS.) | + 
£ -'r,g?K,) AS, e-^ I ^^ - ^ + 1 e-. T (| As, 

Thus, the problem of transport coefficients for specific interaction potential 
(2.2) in our approach is solved. In concluding of this section let us consider 
some limiting cases. 



. 1 Hard spheres potential 



Aeii =0, ASj = 0, 2=1,..., n*, 

(8.8) 
AErj = 0, ASj = 0, j = 1, . . . , m*. 

In this case model MSPI (2.2) transfers to that one for hard spheres, whereas 
kinetic equation (2.13) transfers to that one of the SET (RET) theory. It is 
naturally to expect that results (8.4)-(8.6) should transfer to the well known 
results of the SET theory. This assertion really takes place. With the condition 
(8.8) we have: 
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A 
Hi 
H2 
H2 

X* 



^IgTi^^), 



-A, 

<gTM 



2^o&''(^^) 



Di 
D2 

S(0) 

r(3/2) 



0, 

1, 



TT 



Then 



,{i) 



5(0) 



15„ 



3„cq/ +- 






5n (1 + ^7rna3g^^(ao+)) / m y/^ 



4V2^n^aigT{a^) 



\2kBT 



:8.9) 
^.10) 



p^nkj^T ( 1 + -7rnaj|g-2'^((Tc}"; 



'^^ g^ V ^"^^BTo-;^g'2''(cr+), 



3 5 ( mk-BT' 



1/2 



5 



16 



TT 



^0g2''(f^c[ 



4 x2 

1 + Y^7rna|jg2'^(ff+) 



1/2 



A 



3A;b , 75^A:bTV-^ 1 f^ ^ 2 ^ 

64 V 7rm y 0-0^2(^^0) V 5 



2m 



. 



(8.11) 
(8.12) 

(8.13) 
(8.14) 



Relations (8.12)-(8.14) are obtained from (8.4)-(8.6) with taking into account 
(8.8) are identical to those ones from the SET (RET) theory which are ob- 
tained by means of the standard Chapman- Enskog procedure. 

8.2 Square-well potential 



Ae 



0, 



ASi 



l,...,n*, 



rj 



0, As,- 



0, 



j = 2,...,m* 



AErl ^ 0, 

Asi = f3e. 



^.15) 



In this case initial MSPI (2.2) transfers into "square-well" one of the DRS 
(RDRS) theory. Defining dri = a, Asi =As = {3 ASri = Pe, where e is the 
square- well depth, and taking into account (8.15) one obtains: 



A - ^og2>o) + ^'g2i^-) {e-"^ -1} , 
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i^2 - '^k^i^^) + ^'g?i^') e"^ 5(As), 



//3-^ + e-r(|A.), 
Di^a3g^^(a-)Ase-^^iJ3(As), 

oo 

E;(As) = e^'' — As -2 dx x^Vx^+Ai 






5(0) 



8v^n2 {a* + ^a^gl'^ia-) As^ e-^^} 
5n (l + ^vrnA) + |0FnV3g-^''(a-) Ase^^'HsiAs) / ^ x 1/2 



8v^n2{A* + j^a2g2''(^")As2e-^«} V2A;bT 



and 



p^nA;BT('l + -7rnAV (8.17) 

4 / 

K^-n^^JirmkBTHs, (8.18) 

r/-^^fi: + inA;BT|l + -|v^niJiU(°), (8.19) 

X-^^n+ ^nks (^^] '^' (1 + ^v^nifij a^'\ (8.20) 

2m 4 y m y L 5 J 

Relations for transport coefficients (8.18)-(8.20) coincide with those ones for 
DRS (RDRS) theory which were obtained by means of the standard Chap- 
man-Enskog procedure. Of course, only the first approximation on gradients 
of hydrodynamic parameters is implied everywhere. 



8.3 Smooth long-range potential 

Finally let us consider briefly the case, when 



Aeii ^ 0, Aau ^ 0, i=l,...,n*, n* -^ 00, 
ASrj — > 0, AcTrj — * 0, j = 1, . . . , m*, m* -^ 00, 

and additional condition for (8.21): 
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(8.21) 



_A£k 



A(7rj 



^1 ±1 / 



(8.22) 



From the geometrical point of view, (8.21) and (8.22) correspond to the case 
when MSPI (2.2) is "merged" with some smooth long-range potential 0j at 
r > a. It can be shown that in this case 

p ^ nk^T (^1 + '^nna^gl^a^) - ^vrn^ J dr r3g^^(r)0;(r)) . (8.23) 

Expressions for k, rj and A are completely similar to (8.12)-(8.14) of the SET 
(RET) theory with only difference in the form for g2^{o'o). In SET (RET) 
g2^{o'Q) is some binary equilibrium correlation function of hard spheres on 
contact whereas here it is binary equilibrium correlation function of a system 
with interaction potential of type of hard spheres plus some long "tail" (pt, 
r > o"o- Thus, one obtains final relations for p (8.23) and k, rj and A of the 
KMFT theory [19]. 



9 Numerical calculations 

First of all, let us remember that in the theory under consideration we deal 
with some multistep potential of interaction (2.2). If one has any real infor- 
mation about real (smooth, of course) potential of interaction, it should deal 
with large amount of definition parameters. But really, when interaction po- 
tential is known, the number of independent master parameters is shortened 
greatly. That is the necessary condition, because a model interaction potential 
should approximate real one more or less correctly. First question appearing 
here consists in the representation of an initial smooth interaction potential 
by a multistep one. Let us consider one possible way of definition in which 
all distances between walls of the same kind are equal, i.e.: A an = const, 
i = 1, . . . , n*, AcXrj = const, j = 1, . . . , m*. Then, to define model interaction 
potential one needs (Jq - position of hard sphere wall, cTmax ^ position of the 
most removed attractive wall ((Xmax = ^r„. )) "^p ^ the number of short lengths 
dividing repulsive area [ctq, o-mcan] and nip - the number of short lengths di- 
viding attractive area [o"mcan, o'max], where (Jmcan is the minimum position of 
a real interaction potential. Now MSPI is built. Numbers of repulsive n* and 
attractive m* walls are uniquely determined via numbers of dividing lengths 
Up and rrip. In this representation of some real interaction potential by MSPI, 
one realizes original entwining of model potential around real one. 

Second question is the problem of optimal dividing, i.e. definition of parame- 
ters (To, o"max, np, TUp in such a manner to obtain good results right in the first 
approximation. We tried to solve this problem numerically. 
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Numerical computations of transport coefficients were carried out for Argon 
with the Lennard- Jones potential 



'real 



?LJ 



Ae 



LJ 



^V- f^Y 



r ) 



\ r ) 



(9.1) 



where ctlj = 3.405 A, 5lj/A;b = 119.8 K. 

Starting point in numerical analysis of transport coefficients of our theory are 
relations (8.4)-(8.6) with additional equation for binary equilibrium correla- 
tion function ^^ of a system with potential in a form of multistep function. 
In our calculation we used for ^2 the following approximation: 



i-i^r 



JO) I 



g2'(r)exp{-/30(r)}, 
r) = ip{r), 



(9.2) 



where gg (^) is the binary equilibrium correlation function of hard spheres of 
diameter (Tq. Its analytical expression is well known [27]. 

Table 1 

Parameters for different theories and calculations for transport coefficient rj. Bot- 
tom part contains square displacement of results of SET (RET), MET (BH), DRS 
(RDRS) theories and our theory denoted by GDRS (i.e. generalized DRS) from MD 
simulation. The GDRS result is seen to be the most close to MD simulation. The 
same parameters were used for calculation of another transport coefficients. 



SET (RET) 


SIGMZ0=1.047 


MET (BH) 


1.068 + 0.3837 (A;BT/eLj) 
f^ol j "^LJ 1.000 + 0.4293 (A;BT/eL J ) 


DRS (RDRS) 


SIGMZ0=0.891, SIGMZM=1.342, EZDRS=0.929 


GDRS 


SIGMZ0=0.940, SIGMZM=1.960, np=3, mp=9, n*=2, m*=6 



MD 



SET (RET) MET (BH) DRS (RDRS) GDRS 



0.0 



0.00125 



0.00794 



0.000217 



0.000206 



First, one calculates transport coefficients along gas-liquid saturation curve. 
There were 5 points of calculation {pi = mrii, Tj, i = 1, . . . , 5) along curve of 
saturation for which such transport coefficient as shear viscosity 77 is known 
from the MD simulation [19]. MSPI parameters n^, rup, SIGMZO=cro/crLj, 
SIGMZM=amax/c"o were defined from the minimum of square displacement of 
the theory from corresponding MD results. Parameters of the DRS (RDRS) 
theory were defined in much the same way: SIGMZO=o"o/o"ljj SIGMZM=cr/cro, 
EDRS=£/£lj, as weU as for SET (RET) theory: SIGMZO=ao/aLj. Table 1 
shows the results. Table 2 shows all results of calculation of transport coeffi- 
cients by different theories. Their comparison with experimental data and MD 
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Table 2 

Transport coefficients k, rj and A calculated by different theories. 

Bulk viscosity k, 10^^ Pa sec 

p, g/cm^ SET MET DRS GDRS 

1.4327 0.33387 0.26739 0.43672 0.43371 

1.4180 0.32270 0.25654 0.40946 0.40538 

1.2777 0.22253 0.17126 0.25314 0.24708 

1.1621 0.16222 0.12277 0.17466 0.16928 

0.8017 0.05092 0.03919 0.05914 0.05653 

Shear viscosity rj, 10^^ Pa sec 

p, g/cm3 SET MET DRS GDRS 

0.2970 0.27460 0.22428 0.28794 0.28953 

0.2620 0.26633 0.21627 0.27144 0.27189 

0.1734 0.19113 0.15248 0.17491 0.17248 

0.1255 0.14577 0.11622 0.12627 0.12383 

0.5790 0.06014 0.05210 0.05134 0.05087 

Thermal conductivity A, W/(m K) 

T, K SET MET DRS GDRS 

83.90 0.22107 0.18186 0.17078 0.16850 

86.50 0.21468 0.17566 0.16187 0.15877 

104.50 0.15622 0.12602 0.10790 0.10325 

119.56 0.12080 0.09763 0.08029 0.07603 

147.10 0.05254 0.04608 0.03484 0.03309 

simulations are presented in Figs. 3 and 4. It is clear to see that GDRS results 
practically coinside with experimental data in a wide range of densities and 
temperatures. 

10 Concluding remarks 

Let us discuss areas of application of kinetic equation (2.13). We should re- 
member conditions of general derivation of this equation within the frame of 
Bogolubov-Zubarev approach [29,30]. Specific demand to the geometry of a 
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Fig. 3. Transport coefficients for Argon, a) bulk viscosity k along the liquid- vapour 
curve. X-axis is in units of /9(g/cni^), especially: XAZTI^ 1.4180, 1.1621 and 0.8017 
for 1, 2, 3 and 4, correspondingly, b) shear viscosity t]. x-axis is in units of 
(/)(g/cm3),T(K)), especially: pi = 1.43, Ti = 83.9, p2 = 1.28, Ts = 104.5, ps = 1.16, 
Tg = 119.56 and p^ = 0.802, T4 = 147.1. c) ?? = r]{T) at p = /)„; d) r/ = r]{p) at 
T = 139.7 K. Experimental data plotted in c) and d) are borrowed from [28]. 

potential and to the density of a system is that the mean free path l{ should 
be greatly less than Act - a minimal clearance between the walls. In such a 
way, one should expect that bigger distance between walls and higher density 
give smaller error in the kinetic equation. This error is introduced by the lim- 
iting condition for interaction time |r| -^ +0 [30]. On the other hand, next 
to the theory error there is an error caused by a deviation of the multistep 
potential of interaction from a real one. Real potential is smooth and the error 
is less when clearance between walls is smaller. In the limit (8.21) this error 
is the smallest. One can observe that these two types of errors have oppo- 
site tendencies. So, to apply obtained kinetic equation to systems with real 
smooth interparticle interaction potential in view of a geometry of MSPI one 
should find a compromise solution. Firs of all, MSPI should approximate real 
potential more or less well. At the same time the condition /f <^ Aa must 
obey. This raises the question of optimal dividing of a real potential of inter- 
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Fig. 4. Thermal conductivity A of Argon, a) MD simulations and different theories 
calculations in the same points as in Fig. 3 b). b) A = A(T) at p = 2pcr, different 
theories are compared with experimental data, c) A = A(T) at p = 2pcr', d) A = A(p) 
at T = 298 K. All experimental data plotted in this figure are borrowed from 



action on some multistep one. Density decreasing makes impossible obtaining 
the Boltzmann analogue from the equation under consideration in the limit 
n ^ 0. Let us evaluate numerically. Let ctq is the position of a hard sphere, Aa 
is the minimal distance between walls, fTmax — 2o"o is the location of the most 
removed attractive wall. It is well known from the theory of rarefied gases 
[1,31] that the mean free path /f ?« 1 / \/27ma'^^^. In dense gases it decreases 
in the first approximation in g2*^(crQ") times where g'2^{o'Q) is contact value of 
binary equilibrium correlation function [23]. Thus, l{ ~ l/4:\/27inaQg2^{aQ). 
Introducing dimensionless density A = |7rno"Q, one obtains: 



A(T 



> 



2AV2nAgl\a+) 



7- 



:io.ii 



For A = 0.25 and g2*^(cr(j') ~ 2.5 one obtains 7 ~ 1/25. As far as initial 
preconditions of the theory do not obey then in the limit (8.21), the theory 
error reaches its maximal value. However, kinetic equation transfers then into 
the equation of the kinetic mean field theory [19]. 
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